## Purpose: 

rm(list=ls())
suppressPackageStartupMessages(
  {
    library(tidyverse)
    library(dplyr)
    library(estimatr)
    library(ggplot2)
    library(ggthemes)
    library(texreg)
    library(lfe)
    library(rdrobust)
    library(grid)
    library(gridExtra)
    library(tigris)
    library(sf)
    library(tidycensus)
  })

load("crime_atx_geo.Rdata")
load("stops_atx_geo.Rdata")

#### Creating rdit ####
## Crime
crime_ts_geo = 
  crime_atx_geo_ts %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            violent = mean(violent, na.rm = TRUE))

crime_viol_ts_geo = 
  crime_atx_geo_ts %>% 
  filter(violent == 1) %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE),
            violent = mean(violent, na.rm = TRUE))

crime_ts_geo$blmrv = crime_ts_geo$date - min(crime_ts_geo$date[crime_ts_geo$blm == 1])
crime_viol_ts_geo$blmrv = crime_viol_ts_geo$date - min(crime_viol_ts_geo$date[crime_viol_ts_geo$blm == 1])

rdout_crime_atx = with(crime_ts_geo, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))

## stops
stops_ts_geo = 
  stops_atx_geo_ts %>% 
  mutate(bg = "fips") %>% 
  mutate(count = 1) %>% 
  group_by(bg, date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_ts_geo$blmrv = stops_ts_geo$date - min(stops_ts_geo$date[stops_ts_geo$blm == 1])

rdout_stops_atx = with(stops_ts_geo, rdrobust(x = blmrv, y = count, p = 1, all = TRUE))

##### Separate terciles to calculate total pop within these #####
## crime
crime_atx_geo3 = 
  crime_atx_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high")) %>% 
  mutate(bl_3 = ifelse(bl_low3 == 1, "low", "high")) %>% 
  mutate(bg = as.numeric(fips))

crime_atx_geo_inc3 = 
  crime_atx_geo3 %>%
  group_by(bg) %>% 
  summarise(pop_total = mean(pop_total))

crime_atx_geo_inc3 = 
  crime_atx_geo3 %>%
  group_by(inc_3) %>% 
  summarise(pop_total_inc3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_atx_geo_nw3 = 
  crime_atx_geo3 %>%
  group_by(nw_3) %>%
  summarise(pop_total_nw3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_atx_geo_bl3 = 
  crime_atx_geo3 %>%
  group_by(bl_3) %>%
  summarise(pop_total_bl3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

crime_atx_geo3_ts = crime_atx_geo3
crime_atx_geo_inc3 <- as.data.frame(crime_atx_geo_inc3)
crime_atx_geo_nw3 <- as.data.frame(crime_atx_geo_nw3)
crime_atx_geo_bl3 <- as.data.frame(crime_atx_geo_bl3)
crime_atx_geo3_ts = merge(crime_atx_geo3_ts, crime_atx_geo_inc3, by = "inc_3")
crime_atx_geo3_ts = merge(crime_atx_geo3_ts, crime_atx_geo_nw3, by = "nw_3")
crime_atx_geo3_ts = merge(crime_atx_geo3_ts, crime_atx_geo_bl3, by = "bl_3")

## stops
stops_atx_geo3 = 
  stops_atx_geo_ts %>%
  mutate(inc_3 = ifelse(inc_low3 == 1, "low", "high")) %>%
  mutate(nw_3 = ifelse(nw_low3 == 1, "low", "high")) %>% 
  mutate(bl_3 = ifelse(bl_low3 == 1, "low", "high")) %>% 
  mutate(bg = as.numeric(fips))

stops_atx_geo_inc3 = 
  stops_atx_geo3 %>%
  group_by(bg) %>% 
  summarise(pop_total = mean(pop_total))

stops_atx_geo_inc3 = 
  stops_atx_geo3 %>%
  group_by(inc_3) %>% 
  summarise(pop_total_inc3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stops_atx_geo_nw3 = 
  stops_atx_geo3 %>%
  group_by(nw_3) %>%
  summarise(pop_total_nw3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stops_atx_geo_bl3 = 
  stops_atx_geo3 %>%
  group_by(bl_3) %>%
  summarise(pop_total_bl3 = sum(pop_total, na.rm = TRUE)) %>% 
  mutate(geometry = NULL)

stops_atx_geo3_ts = stops_atx_geo3
stops_atx_geo_inc3 <- as.data.frame(stops_atx_geo_inc3)
stops_atx_geo_nw3 <- as.data.frame(stops_atx_geo_nw3)
stops_atx_geo_bl3 <- as.data.frame(stops_atx_geo_bl3)
stops_atx_geo3_ts = merge(stops_atx_geo3_ts, stops_atx_geo_inc3, by = "inc_3")
stops_atx_geo3_ts = merge(stops_atx_geo3_ts, stops_atx_geo_nw3, by = "nw_3")
stops_atx_geo3_ts = merge(stops_atx_geo3_ts, stops_atx_geo_bl3, by = "bl_3")

#### Create plots ####
## crime
# income
crime_atx_geo_inclow3 = crime_atx_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "A. Crime (Income Lowest Tercile)")

crime_atx_geo_inchigh3 = crime_atx_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "B. Crime (Income Highest Tercile)")

# Nonwhite pop
crime_atx_geo_nwlow3 = crime_atx_geo3_ts %>% 
  filter(nw_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "C. Crime (Nonwhite Population Lowest Tercile)")

crime_atx_geo_nwhigh3 = crime_atx_geo3_ts %>%
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "D. Crime (Nonwhite Population Highest Tercile)")

# Black pop
crime_atx_geo_bllow3 = crime_atx_geo3_ts %>% 
  filter(bl_low3 == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "E. Crime (Black Population Lowest Tercile)")

crime_atx_geo_blhigh3 = crime_atx_geo3_ts %>%
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "F. Crime (Black Population Highest Tercile)")

## Stops
# income
stops_atx_geo_inclow3 = stops_atx_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "A. Stops (Income Lowest Tercile)")

stops_atx_geo_inchigh3 = stops_atx_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "B. Stops (Income Highest Tercile)")

# nonwhite pop
stops_atx_geo_nwlow3 = stops_atx_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "C. Stops (Nonwhite Pop Lowest Tercile)")

stops_atx_geo_nwhigh3 = stops_atx_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "D. Stops (Nonwhite Pop Highest Tercile)")

# black pop
stops_atx_geo_bllow3 = stops_atx_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "E. Stops (Black Pop Lowest Tercile)")

stops_atx_geo_blhigh3 = stops_atx_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) %>% 
  filter(date >= as.Date("2016-01-01")) %>%
  ggplot() + 
  geom_point(aes(x = date, y = count),
             size = .1,
             alpha = .2) +
  geom_smooth(aes(x = date, y = count, group = blm),
              se = FALSE,
              size = .3,
              col = "black") + 
  geom_vline(xintercept = as.Date("2020-05-29"),
             linetype = 2,
             linewidth = .5) + 
  labs(x = "Date", y = "Rate per 100000 people",
       title = "F. Stops (Black Pop Highest Tercile")

## create and save grobs
crime_tercileGrob = arrangeGrob(crime_atx_geo_inclow3, crime_atx_geo_inchigh3, 
                                crime_atx_geo_nwlow3, crime_atx_geo_nwhigh3, 
                                crime_atx_geo_bllow3, crime_atx_geo_blhigh3, 
                                ncol = 3, top=textGrob("Crime by Block Group (Austin)"))

stops_tercileGrob = arrangeGrob(stops_atx_geo_inclow3, stops_atx_geo_inchigh3,
                                stops_atx_geo_nwlow3, stops_atx_geo_nwhigh3,
                                stops_atx_geo_bllow3, stops_atx_geo_blhigh3,
                                ncol = 3, top=textGrob("Stops by Block Group (Austin)"))

ggsave(plot = crime_tercileGrob, filename = "crime_geo_atx.png", width = 16, height = 8)
ggsave(plot = stops_tercileGrob, filename = "stops_geo_atx.png", width = 16, height = 8)

#### Analysis: RDs 25, 50, 100 bw
crime_atx_geo3_nwlow_25 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_low3 == 1)

crime_atx_geo3_nwlow_50 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_low3 == 1)

crime_atx_geo3_nwlow_100 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_low3 == 1)

crime_atx_geo3_nwhigh_25 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_high3 == 1)

crime_atx_geo3_nwhigh_50 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_high3 == 1)

crime_atx_geo3_nwhigh_100 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_high3 == 1)

crime_atx_geo3_inclow_25 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_low3 == 1)

crime_atx_geo3_inclow_50 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_low3 == 1)

crime_atx_geo3_inclow_100 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_low3 == 1)

crime_atx_geo3_inchigh_25 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_high3 == 1)

crime_atx_geo3_inchigh_50 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_high3 == 1)

crime_atx_geo3_inchigh_100 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_high3 == 1)

crime_atx_geo3_bllow_25 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(bl_low3 == 1)

crime_atx_geo3_bllow_50 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(bl_low3 == 1)

crime_atx_geo3_bllow_100 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(bl_low3 == 1)

crime_atx_geo3_blhigh_25 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(bl_high3 == 1)

crime_atx_geo3_blhigh_50 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(bl_high3 == 1)

crime_atx_geo3_blhigh_100 = crime_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(bl_high3 == 1)

stops_atx_geo3_nwlow_25 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_low3 == 1)

stops_atx_geo3_nwlow_50 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_low3 == 1)

stops_atx_geo3_nwlow_100 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_low3 == 1)

stops_atx_geo3_nwhigh_25 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(nw_high3 == 1)

stops_atx_geo3_nwhigh_50 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(nw_high3 == 1)

stops_atx_geo3_nwhigh_100 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(nw_high3 == 1)

stops_atx_geo3_inclow_25 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_low3 == 1)

stops_atx_geo3_inclow_50 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_low3 == 1)

stops_atx_geo3_inclow_100 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_low3 == 1)

stops_atx_geo3_inchigh_25 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(inc_high3 == 1)

stops_atx_geo3_inchigh_50 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(inc_high3 == 1)

stops_atx_geo3_inchigh_100 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(inc_high3 == 1)

stops_atx_geo3_bllow_25 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(bl_low3 == 1)

stops_atx_geo3_bllow_50 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(bl_low3 == 1)

stops_atx_geo3_bllow_100 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(bl_low3 == 1)

stops_atx_geo3_blhigh_25 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-05-04") & date <= as.Date("2020-06-22")) %>% 
  filter(bl_high3 == 1)

stops_atx_geo3_blhigh_50 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-04-08") & date <= as.Date("2020-07-17")) %>% 
  filter(bl_high3 == 1)

stops_atx_geo3_blhigh_100 = stops_atx_geo3_ts %>% 
  filter(date >= as.Date("2020-02-17") & date <= as.Date("2020-09-05")) %>% 
  filter(bl_high3 == 1)

## Separate by outcomes
crime_atx_geo3_inclow_ts = crime_atx_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_atx_geo3_inchigh_ts = crime_atx_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_atx_geo3_nwlow_ts = crime_atx_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_atx_geo3_nwhigh_ts = crime_atx_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_atx_geo3_bllow_ts = crime_atx_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_atx_geo3_blhigh_ts = crime_atx_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_atx_geo3_inclow_ts = stops_atx_geo3_ts %>% 
  filter(inc_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_atx_geo3_inchigh_ts = stops_atx_geo3_ts %>% 
  filter(inc_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_inc3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_atx_geo3_nwlow_ts = stops_atx_geo3_ts %>% 
  filter(nw_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_atx_geo3_nwhigh_ts = stops_atx_geo3_ts %>% 
  filter(nw_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_nw3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_atx_geo3_bllow_ts = stops_atx_geo3_ts %>% 
  filter(bl_low3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_atx_geo3_blhigh_ts = stops_atx_geo3_ts %>% 
  filter(bl_high3 == 1) %>%
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count/(pop_total_bl3/100000), na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

crime_atx_geo3_inclow_ts$blmrv = 
  crime_atx_geo3_inclow_ts$date - min(crime_atx_geo3_inclow_ts$date[crime_atx_geo3_inclow_ts$blm == 1])
crime_atx_geo3_inchigh_ts$blmrv = 
  crime_atx_geo3_inchigh_ts$date - min(crime_atx_geo3_inchigh_ts$date[crime_atx_geo3_inchigh_ts$blm == 1])
crime_atx_geo3_nwlow_ts$blmrv = 
  crime_atx_geo3_nwlow_ts$date - min(crime_atx_geo3_nwlow_ts$date[crime_atx_geo3_nwlow_ts$blm == 1])
crime_atx_geo3_nwhigh_ts$blmrv = 
  crime_atx_geo3_nwhigh_ts$date - min(crime_atx_geo3_nwhigh_ts$date[crime_atx_geo3_nwhigh_ts$blm == 1])
crime_atx_geo3_bllow_ts$blmrv = 
  crime_atx_geo3_bllow_ts$date - min(crime_atx_geo3_bllow_ts$date[crime_atx_geo3_bllow_ts$blm == 1])
crime_atx_geo3_blhigh_ts$blmrv = 
  crime_atx_geo3_blhigh_ts$date - min(crime_atx_geo3_blhigh_ts$date[crime_atx_geo3_blhigh_ts$blm == 1])

stops_atx_geo3_inclow_ts$blmrv = 
  stops_atx_geo3_inclow_ts$date - min(stops_atx_geo3_inclow_ts$date[stops_atx_geo3_inclow_ts$blm == 1])
stops_atx_geo3_inchigh_ts$blmrv = 
  stops_atx_geo3_inchigh_ts$date - min(stops_atx_geo3_inchigh_ts$date[stops_atx_geo3_inchigh_ts$blm == 1])
stops_atx_geo3_nwlow_ts$blmrv = 
  stops_atx_geo3_nwlow_ts$date - min(stops_atx_geo3_nwlow_ts$date[stops_atx_geo3_nwlow_ts$blm == 1])
stops_atx_geo3_nwhigh_ts$blmrv = 
  stops_atx_geo3_nwhigh_ts$date - min(stops_atx_geo3_nwhigh_ts$date[stops_atx_geo3_nwhigh_ts$blm == 1])
stops_atx_geo3_bllow_ts$blmrv = 
  stops_atx_geo3_bllow_ts$date - min(stops_atx_geo3_bllow_ts$date[stops_atx_geo3_bllow_ts$blm == 1])
stops_atx_geo3_blhigh_ts$blmrv = 
  stops_atx_geo3_blhigh_ts$date - min(stops_atx_geo3_blhigh_ts$date[stops_atx_geo3_blhigh_ts$blm == 1])

datasets_atx_geo = list(crime_atx_geo3_inclow_ts %>% as.data.frame %>% mutate(count_crimeatx_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          crime_atx_geo3_inchigh_ts %>% as.data.frame %>% mutate(count_crimeatx_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          crime_atx_geo3_nwlow_ts %>% as.data.frame %>% mutate(count_crimeatx_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          crime_atx_geo3_nwhigh_ts %>% as.data.frame %>% mutate(count_crimeatx_nw_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          crime_atx_geo3_bllow_ts %>% as.data.frame %>% mutate(count_crimeatx_bl_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          crime_atx_geo3_blhigh_ts %>% as.data.frame %>% mutate(count_crimeatx_bl_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          stops_atx_geo3_inclow_ts %>% as.data.frame %>% mutate(count_stopsatx_inc_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          stops_atx_geo3_inchigh_ts %>% as.data.frame %>% mutate(count_stopsatx_inc_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          stops_atx_geo3_nwlow_ts %>% as.data.frame %>% mutate(count_stopsatx_nw_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          stops_atx_geo3_nwhigh_ts %>% as.data.frame %>% mutate(count_stopsatx_nw_high = count) %>% filter(date >= as.Date("2016-01-01")),
                          stops_atx_geo3_bllow_ts %>% as.data.frame %>% mutate(count_stopsatx_bl_low = count) %>% filter(date >= as.Date("2016-01-01")),
                          stops_atx_geo3_blhigh_ts %>% as.data.frame %>% mutate(count_stopsatx_bl_high = count) %>% filter(date >= as.Date("2016-01-01")))

# outcomes

outcomes = c("count_crimeatx_inc_low", "count_crimeatx_inc_high",
             "count_crimeatx_nw_low", "count_crimeatx_nw_high",
             "count_crimeatx_bl_low", "count_crimeatx_bl_high",
             "count_stopsatx_inc_low", "count_stopsatx_inc_high", 
             "count_stopsatx_nw_low", "count_stopsatx_nw_high",
             "count_stopsatx_bl_low", "count_stopsatx_bl_high")

# polynomials

polys = c(0, 1, 2)

# kernels 

kerns = c("tri", "uni", "epa")

# loop for analysis 

dataset_list = as.list(rep(NA, length(datasets_atx_geo)))

for (i in 1:length(datasets_atx_geo)) {
  
  poly_list = as.list(rep(NA, length(polys)))
  
  for (j in 1:length(polys)) {
    
    print(paste0("J Iteration ", j))
    
    kern_list = as.list(rep(NA, length(kerns)))
    
    for (l in 1:length(kerns)) {
      
      print(paste0("L Iteration ", l))
      
      mat_bw = matrix(NA, 100, 7)
      
      for(z in c(25, 50, 100)) {
        
        print(paste0("Z Iteration ", z))
        
        outrd = 
          rdrobust(x = as.numeric(datasets_atx_geo[[i]][, "blmrv"]),
                   y = datasets_atx_geo[[i]][, outcomes[i]],
                   p = polys[j],
                   h = z,
                   kernel = kerns[l])
        
        mat_bw[z, 1] = outrd$coef[3]
        mat_bw[z, 2] = outrd$se[3]
        mat_bw[z, 3] = outrd$pv[3]
        mat_bw[z, 4] = polys[j]
        mat_bw[z, 5] = kerns[l]
        mat_bw[z, 6] = outcomes[i]
        mat_bw[z, 7] = z
        
      }
      
      kern_list[[l]] = mat_bw %>%
        as.data.frame() %>%
        `colnames<-` (c("Coefficient", "SE", "Pvalue", "Polynomial", "Kernel", "Outcome", "BW"))%>% 
        mutate(Coefficient = as.numeric(as.character(Coefficient)),
               SE = as.numeric(as.character(SE)),
               Pvalue = as.numeric(as.character(Pvalue)),
               Polynomial = as.numeric(as.character(Polynomial)),
               Kernel = as.character(Kernel),
               Outcome = as.character(Outcome),
               BW = as.character(BW))
    }
    
    poly_list[[j]] = kern_list
    
  }
  
  dataset_list[[i]] = poly_list
  
}

atx_rdit_geo = dataset_list %>% 
  unlist(recursive = FALSE) %>%
  unlist(recursive = FALSE) %>%
  do.call(rbind.data.frame, .) %>% 
  mutate(Outcome = factor(Outcome,
                          levels = c("count_crimeatx_inc_low", "count_crimeatx_inc_high",
                                     "count_crimeatx_nw_low", "count_crimeatx_nw_high",
                                     "count_crimeatx_bl_low", "count_crimeatx_bl_high",
                                     "count_stopsatx_inc_low", "count_stopsatx_inc_high", 
                                     "count_stopsatx_nw_low", "count_stopsatx_nw_high",
                                     "count_stopsatx_bl_low", "count_stopsatx_bl_high"),
                          labels = c("A. Crime, Low Income",
                                     "B. Crime, High Income",
                                     "C. Crime, Low Percentage Nonwhite",
                                     "D. Crime, High Percentage Nonwhite",
                                     "E. Crime, Low Percentage Black",
                                     "F. Crime, High Percentage Black",
                                     "G. Terry Stops, Low Income",
                                     "H. Terry Stops, High Income",
                                     "I. Terry Stops, Low Percentage Nonwhite",
                                     "J. Terry Stops, High Percentage Nonwhite",
                                     "K. Terry Stops, Low Percentage Black",
                                     "L. Terry Stops, High Percentage Black"
                          )),
         Polynomial = factor(Polynomial, levels = c(0, 1, 2),
                             labels = c("DIM", "Linear", "Quadratic")),
         Kernel = factor(Kernel, levels = c("tri", "uni", "epa"),
                         labels = c("Triangular",
                                    "Uniform",
                                    "Epanechnikov"))) %>% 
  filter(!is.na(BW))

atx_rdit_geo = atx_rdit_geo %>% 
  mutate(kern_poly = factor(paste0(Kernel, ", ", Polynomial), levels = c("Triangular, DIM", 
                                                                         "Triangular, Linear",
                                                                         "Triangular, Quadratic",
                                                                         "Uniform, DIM",
                                                                         "Uniform, Linear",
                                                                         "Uniform, Quadratic",
                                                                         "Epanechnikov, DIM",
                                                                         "Epanechnikov, Linear",
                                                                         "Epanechnikov, Quadratic")))

atx_rdit_geo$BW = as.numeric(atx_rdit_geo$BW)

#### Coefficient tests of difference ####
## Create function for s
s <- function(n1, n2, p, s1, s2){sqrt(((n1 - p)*((s1)^2) + (n2- p)*((s2)^2)/(n1 + n2 - 2*p)))}

se <- function(s, SEb1, SEb2, s1, s2){s*sqrt(((SEb1/s1)^2) + (SEb2/s2)^2)}

z <- function(b1, b2, se){(b1-b2)/se}

b_est <- function(b1, b2){(b1-b2)}

s_crimeatx_25_inc <- s(n1 = 3706, n2 = 4924, p = 2, s1 = 7.746181e-04, s2 = 3.130744e-04)
se_catx_25_inc <- se(s = s_crimeatx_25_inc, SEb1 = 7.746181e-04, SEb2 = 3.130744e-04, s1 = 7.746181e-04, s2 = 3.130744e-04)
z_catx_25_inc <- z(b1 = -1.102393e-04, b2 = -2.143327e-04, se = se_catx_25_inc)
b_catx_25_inc <- b_est(b1 = -1.102393e-04, b2 = -2.143327e-04)
p_catx_25_inc <- 2*pnorm(abs(z_catx_25_inc), mean = 0, lower.tail = FALSE)

s_crimeatx_25_nw <- s(n1 = 2612, n2 = 5324, p = 2, s1 = 9.679808e-04, s2 = 2.288060e-04)
se_catx_25_nw <- se(s = s_crimeatx_25_nw, SEb1 = 9.679808e-04, SEb2 = 2.288060e-04, s1 = 9.679808e-04, s2 = 2.288060e-04)
z_catx_25_nw <- z(b1 = -1.081278e-04, b2 = -1.125114e-04, se = se_catx_25_nw)
b_catx_25_nw <- b_est(b1 = -1.081278e-04, b2 = -1.125114e-04)
p_catx_25_nw <- 2*pnorm(abs(z_catx_25_nw), mean = 0, lower.tail = FALSE)

s_crimeatx_25_bl <- s(n1 = 2748, n2 = 4702, p = 2, s1 = 7.867861e-04, s2 = 1.830940e-04)
se_catx_25_bl <- se(s = s_crimeatx_25_bl, SEb1 = 7.867861e-04, SEb2 = 1.830940e-04, s1 = 7.867861e-04, s2 = 1.830940e-04)
z_catx_25_bl <- z(b1 = -4.342413e-04, b2 = 3.673437e-04, se = se_catx_25_bl)
b_catx_25_bl <- b_est(b1 = -4.342413e-04, b2 = 3.673437e-04)
p_catx_25_bl <- 2*pnorm(abs(z_catx_25_bl), mean = 0, lower.tail = FALSE)

s_crimeatx_50_inc <- s(n1 = 7270, n2 = 9730, p = 2, s1 = 5.049419e-04, s2 = 2.253926e-04)
se_catx_50_inc <- se(s = s_crimeatx_50_inc, SEb1 = 5.049419e-04, SEb2 = 2.253926e-04, s1 = 5.049419e-04, s2 = 2.253926e-04)
z_catx_50_inc <- z(b1 = -3.361063e-04, b2 = -3.409916e-04, se = se_catx_50_inc)
b_catx_50_inc <- b_est(b1 = -3.361063e-04, b2 = -3.409916e-04)
p_catx_50_inc <- 2*pnorm(abs(z_catx_50_inc), mean = 0, lower.tail = FALSE)

s_crimeatx_50_nw <- s(n1 = 5085, n2 = 10957, p = 2, s1 = 6.338638e-04, s2 = 1.629000e-04)
se_catx_50_nw <- se(s = s_crimeatx_50_nw, SEb1 = 6.338638e-04, SEb2 = 1.629000e-04, s1 = 6.338638e-04, s2 = 1.629000e-04)
z_catx_50_nw <- z(b1 = -3.924567e-04, b2 = -2.353951e-04, se = se_catx_50_nw)
b_catx_50_nw <- b_est(b1 = -3.924567e-04, b2 = -2.353951e-04)
p_catx_50_nw <- 2*pnorm(abs(z_catx_50_nw), mean = 0, lower.tail = FALSE)

s_crimeatx_50_bl <- s(n1 = 5428, n2 = 9526, p = 2, s1 = 5.981170e-04, s2 = 1.308153e-04)
se_catx_50_bl <- se(s = s_crimeatx_50_bl, SEb1 = 5.981170e-04, SEb2 = 1.308153e-04, s1 = 5.981170e-04, s2 = 1.308153e-04)
z_catx_50_bl <- z(b1 = -7.869802e-04, b2 = 1.899914e-04, se = se_catx_50_bl)
b_catx_50_bl <- b_est(b1 = -7.869802e-04, b2 = 1.899914e-04)
p_catx_50_bl <- 2*pnorm(abs(z_catx_50_bl), mean = 0, lower.tail = FALSE)

s_crimeatx_100_inc <- s(n1 = 14859, n2 = 19611, p = 2, s1 = 3.203085e-04, s2 = 1.564996e-04)
se_catx_100_inc <- se(s = s_crimeatx_100_inc, SEb1 = 3.203085e-04, SEb2 = 1.564996e-04, s1 = 3.203085e-04, s2 = 1.564996e-04)
z_catx_100_inc <- z(b1 = -5.314628e-05, b2 = -1.853590e-04, se = se_catx_100_inc)
b_catx_100_inc <- b_est(b1 = -5.314628e-05, b2 = -1.853590e-04)
p_catx_100_inc <- 2*pnorm(abs(z_catx_100_inc), mean = 0, lower.tail = FALSE)

s_crimeatx_100_nw <- s(n1 = 10489, n2 = 21737, p = 2, s1 = 4.143981e-04, s2 = 1.153775e-04)
se_catx_100_nw <- se(s = s_crimeatx_100_nw, SEb1 = 4.143981e-04, SEb2 = 1.153775e-04, s1 = 4.143981e-04, s2 = 1.153775e-04)
z_catx_100_nw <- z(b1 = -4.776321e-04, b2 = -1.658988e-04, se = se_catx_100_nw)
b_catx_100_nw <- b_est(b1 = -4.776321e-04, b2 = -1.658988e-04)
p_catx_100_nw <- 2*pnorm(abs(z_catx_100_nw), mean = 0, lower.tail = FALSE)

s_crimeatx_100_bl <- s(n1 = 11275, n2 = 19194, p = 2, s1 = 4.339007e-04, s2 = 9.156318e-05)
se_catx_100_bl <- se(s = s_crimeatx_100_bl, SEb1 = 4.339007e-04, SEb2 = 9.156318e-05, s1 = 4.339007e-04, s2 = 9.156318e-05)
z_catx_100_bl <- z(b1 = -5.383571e-04, b2 = 5.710699e-05, se = se_catx_100_bl)
b_catx_100_bl <- b_est(b1 = -5.383571e-04, b2 = 5.710699e-05)
p_catx_100_bl <- 2*pnorm(abs(z_catx_100_bl), mean = 0, lower.tail = FALSE)

s_stopsatx_25_inc <- s(n1 = 784, n2 = 1805, p = 2, s1 = 8.979375e-02, s2 = 6.498886e-02)
se_satx_25_inc <- se(s = s_stopsatx_25_inc, SEb1 = 8.979375e-02, SEb2 = 6.498886e-02, s1 = 8.979375e-02, s2 = 6.498886e-02)
z_satx_25_inc <- z(b1 = -4.269214e-01, b2 = -2.943653e-01, se = se_satx_25_inc)
b_satx_25_inc <- b_est(b1 = -4.269214e-01, b2 = -2.943653e-01)
p_satx_25_inc <- 2*pnorm(abs(z_satx_25_inc), mean = 0, lower.tail = FALSE)

s_stopsatx_25_nw <- s(n1 = 637, n2 = 739, p = 2, s1 = 7.495753e-02, s2 = 2.664070e-02)
se_satx_25_nw <- se(s = s_stopsatx_25_nw, SEb1 = 7.495753e-02, SEb2 = 2.664070e-02, s1 = 7.495753e-02, s2 = 2.664070e-02)
z_satx_25_nw <- z(b1 = -2.797799e-01, b2 = -9.744848e-02, se = se_satx_25_nw)
b_satx_25_nw <- b_est(b1 = -2.797799e-01, b2 = -9.744848e-02)
p_satx_25_nw <- 2*pnorm(abs(z_satx_25_nw), mean = 0, lower.tail = FALSE)

s_stopsatx_25_bl <- s(n1 = 530, n2 = 1128, p = 2, s1 = 6.569452e-02, s2 = 4.316895e-02)
se_satx_25_bl <- se(s = s_stopsatx_25_bl, SEb1 = 6.569452e-02, SEb2 = 4.316895e-02, s1 = 6.569452e-02, s2 = 4.316895e-02)
z_satx_25_bl <- z(b1 = -5.526287e-01, b2 = -1.021991e-01, se = se_satx_25_bl)
b_satx_25_bl <- b_est(b1 = -5.526287e-01, b2 = -1.021991e-01)
p_satx_25_bl <- 2*pnorm(abs(z_satx_25_bl), mean = 0, lower.tail = FALSE)

s_stopsatx_50_inc <- s(n1 = 1291, n2 = 2763, p = 2, s1 = 6.281433e-02, s2 = 3.872724e-02)
se_satx_50_inc <- se(s = s_stopsatx_50_inc, SEb1 = 6.281433e-02, SEb2 = 3.872724e-02, s1 = 6.281433e-02, s2 = 3.872724e-02)
z_satx_50_inc <- z(b1 = -4.056551e-01, b2 = -2.889409e-01, se = se_satx_50_inc)
b_satx_50_inc <- b_est(b1 = -4.056551e-01, b2 = -2.889409e-01)
p_satx_50_inc <- 2*pnorm(abs(z_satx_50_inc), mean = 0, lower.tail = FALSE)

s_stopsatx_50_nw <- s(n1 = 1109, n2 = 1161, p = 2, s1 = 5.084320e-02, s2 = 1.573426e-02)
se_satx_50_nw <- se(s = s_stopsatx_50_nw, SEb1 = 5.084320e-02, SEb2 = 1.573426e-02, s1 = 5.084320e-02, s2 = 1.573426e-02)
z_satx_50_nw <- z(b1 = -2.915648e-01, b2 = -9.524029e-02, se = se_satx_50_nw)
b_satx_50_nw <- b_est(b1 = -2.915648e-01, b2 = -9.524029e-02)
p_satx_50_nw <- 2*pnorm(abs(z_satx_50_nw), mean = 0, lower.tail = FALSE)

s_stopsatx_50_bl <- s(n1 = 805, n2 = 1858, p = 2, s1 = 4.554690e-02, s2 = 2.603792e-02)
se_satx_50_bl <- se(s = s_stopsatx_50_bl, SEb1 = 4.554690e-02, SEb2 = 2.603792e-02, s1 = 4.554690e-02, s2 = 2.603792e-02)
z_satx_50_bl <- z(b1 = -4.889987e-01, b2 = -1.289459e-01, se = se_satx_50_bl)
b_satx_50_bl <- b_est(b1 = -4.889987e-01, b2 = -1.289459e-01)
p_satx_50_bl <- 2*pnorm(abs(z_satx_50_bl), mean = 0, lower.tail = FALSE)

s_stopsatx_100_inc <- s(n1 = 2790, n2 = 5584, p = 2, s1 = 4.078780e-02, s2 = 2.356830e-02)
se_satx_100_inc <- se(s = s_stopsatx_100_inc, SEb1 = 4.078780e-02, SEb2 = 2.356830e-02, s1 = 4.078780e-02, s2 = 2.356830e-02)
z_satx_100_inc <- z(b1 = -2.871124e-01, b2 = -1.892004e-01, se = se_satx_100_inc)
b_satx_100_inc <- b_est(b1 = -2.871124e-01, b2 = -1.892004e-01)
p_satx_100_inc <- 2*pnorm(abs(z_satx_100_inc), mean = 0, lower.tail = FALSE)

s_stopsatx_100_nw <- s(n1 = 2711, n2 = 2904, p = 2, s1 = 3.322483e-02, s2 = 9.418335e-03)
se_satx_100_nw <- se(s = s_stopsatx_100_nw, SEb1 = 3.322483e-02, SEb2 = 9.418335e-03, s1 = 3.322483e-02, s2 = 9.418335e-03)
z_satx_100_nw <- z(b1 = -2.066181e-01, b2 = -5.697993e-02, se = se_satx_100_nw)
b_satx_100_nw <- b_est(b1 = -2.066181e-01, b2 = -5.697993e-02)
p_satx_100_nw <- 2*pnorm(abs(z_satx_100_nw), mean = 0, lower.tail = FALSE)

s_stopsatx_100_bl <- s(n1 = 1862, n2 = 3584, p = 2, s1 = 3.075177e-02, s2 = 1.602009e-02)
se_satx_100_bl <- se(s = s_stopsatx_100_bl, SEb1 = 3.075177e-02, SEb2 = 1.602009e-02, s1 = 3.075177e-02, s2 = 1.602009e-02)
z_satx_100_bl <- z(b1 = -2.603917e-01, b2 = -9.695119e-02, se = se_satx_100_bl)
b_satx_100_bl <- b_est(b1 = -2.603917e-01, b2 = -9.695119e-02)
p_satx_100_bl <- 2*pnorm(abs(z_satx_100_bl), mean = 0, lower.tail = FALSE)
